% This function computes market share derivatives w.r.t. prices

function [D_s1_p1, D_s2_p1, D_s2_p2, D_s1_ind_p1]...
    = D_s_price(p1, s1_predicted_y_c, s1_predicted_y_m, pr_s, ...
    resource_y1, resource_y2, pr_y, N_y, ...
    NN, K, M,...
    myopic, delta, B1_c, B2_c, beta_c, T1, alpha, rho, crra, sigma_c1)

KM = K*M;

%---------------------
% D_s1_p1: NN x 1. (j)th element = derivative of s_j1  w.r.t. p_j1

% 1st-period consumption after paying premiums, N_y x NN
c1 = repmat(resource_y1,1,NN) - repmat(p1',N_y,1);

% 2nd-period consumption when not lapse & the price is set to p1 
c2 = repmat(resource_y2,1,NN) - repmat(p1',N_y,1);

% for each insurer, compute the expected probability of consumers staying.
pr_stay_fixed1 = (1-unique(delta))*ones(KM,1); 
agg_stay = pr_s*pr_stay_fixed1; % NN x 1

% D_v_p1: N_y x NN 
D_v_p1_c = -B1_c*alpha*D_u_c(c1, rho, crra)/sigma_c1;
D_v_p1_m = (-B1_c*alpha*D_u_c(c1, rho, crra)...
    - ((beta_c^T1)*B2_c*alpha)*D_u_c(c2, rho, crra).*repmat(agg_stay', N_y,1))/sigma_c1;

% D_s1_ind_p1: N_y x NN
D_s1_ind_p1_c = D_v_p1_c.*s1_predicted_y_c.*(1-s1_predicted_y_c);
D_s1_ind_p1_m = D_v_p1_m.*s1_predicted_y_m.*(1-s1_predicted_y_m);
D_s1_ind_p1 = myopic*D_s1_ind_p1_m + (1-myopic)*D_s1_ind_p1_c;

% D_s1_p1: NN x 1.
D_s1_p1 = sum(D_s1_ind_p1.*repmat(pr_y, 1, NN))';

%---------------------
% D_s2_p1: NN x K. (jk)th element = derivative of s_jk2 w.r.t. p_j1.
pr_stay_fixed2 = (1-unique(delta))*ones(N_y,KM); 

D_s2_p1 = zeros(NN, KM);

for jj=1:NN
    
    % D_s1_ind_p1_jj: N_y x 1
    D_s1_ind_p1_jj = D_s1_ind_p1(:,jj);
    
    % 1 x K
    D_s2_p1(jj,:) = sum(repmat(D_s1_ind_p1_jj,1,KM).*pr_stay_fixed2.*repmat(pr_y, 1, KM));
   
end

%---------------------
% D_s2_p2: NN x K. (jk)th element = derivative of s_jk2 w.r.t. p_jk2 ->
% zero with exogenous lapses
D_s2_p2 = zeros(NN, KM);
